from pylab import *

######################################################################
#
# TO BE TESTED AGAINST THE ORIGINAL MATLAB VERSION
#
######################################################################

#
# RFCFILTER - Rainflow filter a signal
# This code is based on the rfcfilter.m file contained in the
# WAFO toolbox, created by P Johannesson in 1999, and released under
# the GNU GPL.
#

def rfcfilter(x, h, defi = 0):
  """rfcfilter()
  Filter a time-based signal using the rainflow counting technique.
  The resulting data can then be used to construct a Markov matrix
  that can be used for damage (life) calculations.
  """
  y = zeros((x.shape), dtype = float64)
  m = x.shape[0]
  n = x.shape[1]
  
  y[0,:] = x[0,:]
  j = 0
  if n == 1:
    t0 = 0
    y0 = y[0]
  else:
    t0 = y[0,0]
    y0 = y[0,1]
  
  z0 = 0
  
  # Start of the rainflow filter
  for i in range(1, m):
    fpi = y0 + h
    fmi = y0 - h
    if n == 1:
      ti = 0
      xi = x[i]
    else:
      ti = x[i,0]
      xi = x[i,1]
    
    if z0 == 0:
      if defi == 0:
        test1 = (xi <= fmi)
        test2 = (xi >= fpi)
      else:
        test1 = (xi < fmi)
        test2 = (xi > fpi)
      if test1:
        z1 = -1
      elif test2:
        z1 = +1
      else:
        z1 = 0
      if z1 == 0:
        t1 = t0
        y1 = y0
      else:
        t1 = ti
        y1 = xi
    else:
      if defi == 0:
        test1 = (((z0==1) and (xi<=fmi)) or ((z0==-1) and (xi<fpi)))
        test2 = (((z0==1) and (xi>fmi)) or ((z0==-1) and (xi>=fpi)))
      else:
        test1 = (((z0==+1) and (xi<fmi))  or ((z0==-1) and (xi<=fpi)))
        test2 = (((z0==+1) and (xi>=fmi)) or ((z0==-1) and (xi>fpi)))
      
      #Update z1
      if test1:
        z1 = -1
      elif test2:
        z1 = +1
      else:
        print 'Something is wrong, i = %d' % (i,)
        
      #Update y1
      if z1 != z0:
        t1 = ti
        y1 = xi
      elif z1 == -1:
        #y1 = min([y0 xi])
        if y0 < xi:
          t1 = t0
          y1 = y0
        else:
          t1 = ti
          y1 = xi
      elif z1 == +1:
        #y1 = max([y0 xi])
        if y0 > xi:
          t1 = t0
          y1 = y0
        else:
          t1 = ti
          y1 = xi
    #Update y if y0 is a turning point
    if abs(z0-z1) == 2:
      j = j + 1
      if n == 1:
        y[j] = y0
      else:
        y[j,:] = [t0, y0]
  
    #Update t0, y0, z0
    t0 = t1
    y0 = y1
    z0 = z1
  
  #Update y if last y0 is greater than (or equal) threshold
  if defi == 0:
    test = (abs(y0-y[j]) > h)
  else: #def == 1
    test = (abs(y0-y[j]) >= h)

  if test:
    j=j+1
    if n == 1:
      y[j] = y0
    else:
      y[j,:] = [t0, y0]
  #Truncate y
  y=y[1:j,:]
  return y


from scipy.signal import butter, lfilter
# Calculate Butterworth filter coefficients.
# This will be used to low-pass filter a random signal
[B,A] = butter(8, 0.5)
seed(3)
x = randn(50000,1) * 800.
x1 = lfilter(B, A, x, axis = 0)
y = rfcfilter(x,100.)

